Comparison of the snia_toy_06 model

In [139]:
import os, itertools
import pandas as pd
import numpy as np

from bokeh.plotting import figure, output_file, show
from bokeh.io import output_notebook
from bokeh.palettes import Dark2_5, colorblind
from glob import glob

from bokeh.layouts import layout, row
from bokeh.models.widgets import Tabs, Panel
from bokeh.io import curdoc
from bokeh.plotting import figure
from scipy import interpolate
In [148]:
def read_magnitudes(fname):
    COL_NAMES = ['time', 'U', 'B', 'V', 'R', 'I', 'J', 'H', 'K']
    mag_table = pd.read_csv(fname, delim_whitespace=True, comment='#', names=COL_NAMES)
    mag_table.iloc[:, 1:] = mag_table.iloc[:,1:].replace(99.999, np.nan)
    return mag_table
    
def read_spectrum(fname):
    col_names = ['wavelength']
    with open(fname, encoding='utf8') as fh:
        for i, line in enumerate(fh):
            if i==2:
                times = list(map(np.float64, line.strip().split(':')[1].split()))
    
    spectrum = pd.read_csv(fname, delim_whitespace=True, comment='#', header=None, names=col_names + times)
    return spectrum

def read_tgas(fname):
    col_names = ['velocity']
    with open(fname, encoding='utf8') as fh:
        for i, line in enumerate(fh):
            if i==2:
                times = list(map(np.float64, line.strip().split(':')[1].split()))
    
    spectrum = pd.read_csv(fname, delim_whitespace=True, comment='#', header=None, names=col_names + times)
    return spectrum
In [88]:
# uncomment the following line to download the code comparison stuff
#!wget http://opensupernova.org/~wkerzend/files/ccompare1_models.tgz
#!tar zxf ccompare1_models.tgz
In [85]:
CODE_NAMES = np.sort([os.path.basename(fname).split('_')[-1][:-4] for fname in glob('toy06/lbol_edep_toy06_*.txt')])
CODE_COLORS = {name:color for name, color in zip(CODE_NAMES, colorblind['Colorblind'][8])}

Comparison of pseudo bolometric luminosity

Integration of the synthetic spectra excluding gamma rays ($>50 Angstrom$)

In [147]:
output_notebook()


# create a new plot with a title and axis labels
fig_lbol  = figure(title="Comparison of Lbol", x_axis_label='time [days]', y_axis_label='Log Lbol')
fig_edep  = figure(title="Comparison of energy deposition", x_axis_label='time [days]', y_axis_label='Log Edeptot')

log_lbol_data = {}
for fname in glob('toy06/lbol_edep_toy06_*.txt'):
    name = os.path.basename(fname).split('_')[-1][:-4]
    lbol = pd.read_csv(fname, delim_whitespace=True, comment='#')
    lbol.columns = ['time', 'lbol', 'edep']
    lbol['log_lbol'] = np.log10(lbol.lbol)
    lbol['log_edep'] = np.log10(lbol.edep)
    #lbol.plot_bokeh()
    # add a line renderer with legend and line thickness
    fig_lbol.line(lbol.time, lbol.log_lbol, legend=name, line_width=2, color=CODE_COLORS[name])
    fig_edep.line(lbol.time, lbol.log_edep, legend=name, line_width=2, color=CODE_COLORS[name])
    log_lbol_data[name] = lbol

fig_lbol.legend.location = "top_right"
fig_lbol.legend.click_policy="hide"
fig_edep.legend.location = "top_right"
fig_edep.legend.click_policy="hide"

# show the results
show(row([fig_lbol, fig_edep]))#
Loading BokehJS ...
/Users/wkerzend/miniconda/envs/sn-rad-trans/lib/python3.7/site-packages/pandas/core/series.py:853: RuntimeWarning: divide by zero encountered in log10
  result = getattr(ufunc, method)(*inputs, **kwargs)

Comparing gas temperatures

In [155]:
## Reading the data
tgas_data = {}
for fname in glob('toy06/tgas_toy06_*'):
    name = os.path.basename(fname).split('_')[-1][:-4]
    try:
        tgas_data[name] = read_tgas(fname)
    except UnicodeDecodeError:
        pass
    except IndexError:
        pass
In [158]:
TGAS_TIMES = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400]
output_notebook()
interpolators = {}
for name in CODE_NAMES:
    if name not in tgas_data:
        continue
    interpolators[name] = interpolate.interp1d(tgas_data[name].columns[1:], 
                                               tgas_data[name].iloc[:, 1:].values, 
                                               fill_value=np.nan, bounds_error=False)
all_tabs = []
for time in TGAS_TIMES:
    fig = figure(x_axis_label='velocity [km/s]', y_axis_label='Tgas [K]')
    for name, interpolator in interpolators.items():
        velocity = tgas_data[name].velocity
        fig.line(velocity, interpolator(time), legend=name, line_width=2, color=CODE_COLORS[name])
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title='t={0:.0f}'.format(time))
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...

Comparing the magnitudes

In [92]:
## Reading the data
mag_data = {}
for fname, color in zip(glob('toy06/mags_toy06_*'), colors):
    name = os.path.basename(fname).split('_')[-1][:-4]
    mag_data[name] = read_magnitudes(fname)
In [93]:
output_notebook()

MAG_NAMES = list(mag_data.values())[0].columns[1:]
all_tabs = []
for mag in MAG_NAMES:
    fig = figure()
    for name, mag_table in mag_data.items():
        fig.line(mag_table.time, mag_table[mag], legend=name, line_width=2, color=CODE_COLORS[name])
        fig.y_range.flipped = True
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title=mag)
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...

Comparing the colors

In [95]:
output_notebook()

MAG_NAMES = list(mag_data.values())[0].columns[1:]
all_tabs = []
for mag1, mag2 in zip(MAG_NAMES[:-1], MAG_NAMES[1:]):
    fig = figure()
    for name, mag_table in mag_data.items():
        fig.line(mag_table.time, mag_table[mag1] - mag_table[mag2], legend=name, line_width=2, color=CODE_COLORS[name])
        fig.y_range.flipped = True
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title='{0}-{1}'.format(mag1, mag2))
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...

Comparing the spectra

In [117]:
## Reading the data
spectral_data = {}
for fname, color in zip(glob('toy06/spectra_toy06_*'), colors):
    name = os.path.basename(fname).split('_')[-1][:-4]
    try:
        spectral_data[name] = read_spectrum(fname)
    except UnicodeDecodeError:
        pass
    except IndexError:
        pass
In [141]:
interpolate.interp1d?
In [143]:
SPECTRAL_TIMES = [5, 10, 15, 20, 30, 50, 100, 200, 300, 400]
output_notebook()
interpolators = {}
for name in CODE_NAMES:
    if name not in spectral_data:
        continue
    interpolators[name] = interpolate.interp1d(spectral_data[name].columns[1:], 
                                               spectral_data[name].iloc[:, 1:].values, 
                                               fill_value=np.nan, bounds_error=False)

all_tabs = []
for time in SPECTRAL_TIMES:
    fig = figure()
    for name, interpolator in interpolators.items():
        wavelength = spectral_data[name].wavelength
        fig.line(wavelength, interpolator(time), legend=name, line_width=2, color=CODE_COLORS[name])
    fig.legend.location = "top_right"
    fig.legend.click_policy="hide"
    lout = layout([[fig]])
    tab = Panel(child=lout,title='t={0:.0f}'.format(time))
    all_tabs.append(tab)
        
show(Tabs(tabs=all_tabs))
Loading BokehJS ...
In [ ]: